##########################################################################################
# Title: Figure 3b — (Pop 10–20 km) / (Pop within 10 km) of a BD, by City & Year
# Section: 2 (Figures)
# Purpose: Benchmark relative concentration far vs. near BDs across 25 cities
# IO:
#   IN : Processed_Data/Intermediate/Section_2/Processed_BD_Analysis/
#        ├─ ALL_CITIES_100K_09_2000.rds
#        ├─ ALL_CITIES_100K_09_2010.rds
#        └─ ALL_CITIES_100K_09_2020.rds
#   OUT: Graphs/Figure_3b.png
##########################################################################################

rm(list = ls()); options(stringsAsFactors = FALSE)

suppressPackageStartupMessages({
  library(data.table)
  library(ggplot2)
  library(grid)     # unit()
})

# ---- Project root (align with run_all.R) -------------------------------------
if (!requireNamespace("here", quietly = TRUE)) install.packages("here", quiet = TRUE)
ROOT     <- normalizePath(here::here(), winslash = "/")
PDIR     <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2")
IN_DIR   <- file.path(PDIR, "Processed_BD_Analysis")
GRAPH_DIR<- file.path(ROOT, "Graphs")
dir.create(GRAPH_DIR, recursive = TRUE, showWarnings = FALSE)

# ---- Inputs ------------------------------------------------------------------
f2000 <- file.path(IN_DIR, "ALL_CITIES_100K_09_2000.rds")
f2010 <- file.path(IN_DIR, "ALL_CITIES_100K_09_2010.rds")
f2020 <- file.path(IN_DIR, "ALL_CITIES_100K_09_2020.rds")
stopifnot(file.exists(f2000), file.exists(f2010), file.exists(f2020))

out_file <- file.path(GRAPH_DIR, "Figure_3b.png")

# ---- Helper: compute ratio (15–20 km) / (0–10 km) by city --------------------
compute_ratios <- function(path, year_label) {
  dt <- readRDS(path) |> as.data.table()
  # keep only needed bands and clean
  dt_sub <- dt[bounds %in% c(5L,10L,15L,20L) & !is.na(city_name) & !is.na(pop_sum)]
  # total pop per city x band
  bands  <- dt_sub[, .(pop_total = sum(pop_sum, na.rm = TRUE)), by = .(city_name, bounds)]
  # ensure all bands exist for every city
  all_lvls <- CJ(city_name = unique(bands$city_name), bounds = c(5L,10L,15L,20L))
  bands <- merge(all_lvls, bands, by = c("city_name","bounds"), all.x = TRUE)
  bands[is.na(pop_total), pop_total := 0]
  # wide + ratio
  wide <- dcast(bands, city_name ~ bounds, value.var = "pop_total")
  setnames(wide, c("5","10","15","20"), c("b5","b10","b15","b20"))
  wide[, ratio_10_20 := fifelse(b5 + b10 > 0, (b15 + b20) / (b5 + b10), NA_real_)]
  wide[, year := year_label]
  wide[, .(city_name, ratio_10_20, year)]
}

# ---- Load & stack ------------------------------------------------------------
d2000 <- compute_ratios(f2000, "2000")
d2010 <- compute_ratios(f2010, "2010")
d2020 <- compute_ratios(f2020, "2020")
df_all <- rbindlist(list(d2000, d2010, d2020), use.names = TRUE)

# ---- Prep for plotting -------------------------------------------------------
plot_df <- df_all[city_name != "Mumbai" & !is.na(ratio_10_20)]

# Order cities by 2020 (desc), then reverse factor so largest at top
order_vec <- plot_df[year == "2020"][order(-ratio_10_20), unique(city_name)]
plot_df[, city_name := factor(city_name, levels = rev(order_vec))]
plot_df[, year := factor(year, levels = c("2000","2010","2020"))]

# ---- Colors (Stata-like navy ramp) ------------------------------------------
color_map <- c("2000"="#B0C4DE", "2010"="#4682B4", "2020"="#0B3C6F")

# ---- Axis range --------------------------------------------------------------
ymax  <- max(plot_df$ratio_10_20, na.rm = TRUE)
upper <- if (ymax <= 1.5) 1.5 else ceiling(ymax * 20) / 20  # step ~0.05 above max

# ---- Plot --------------------------------------------------------------------
p <- ggplot(plot_df, aes(x = city_name, y = ratio_10_20)) +
  geom_col(aes(fill = year), width = 0.70, position = position_dodge(width = 0.80)) +
  coord_flip() +
  scale_fill_manual(
    values = color_map, name = NULL,
    guide = guide_legend(
      direction = "horizontal",
      nrow = 1,
      label.position = "bottom",
      keywidth = 3.0, keyheight = 0.35, label.hjust = 0.5
    )
  ) +
  labs(
    y = "(Pop. 10–20 km) / (Pop. within 10 km) of a Business District",
    x = NULL
  ) +
  scale_y_continuous(
    limits = c(0, upper),
    breaks = seq(0, max(1.5, upper), by = 0.5),
    expand = expansion(mult = c(0, 0.05))
  ) +
  geom_hline(yintercept = c(0.5, 1.0), linetype = "dotted", color = "grey60") +
  theme_classic(base_size = 10) +
  theme(
    legend.position      = "bottom",
    legend.justification = "center",
    legend.box           = "vertical",
    legend.box.spacing   = unit(0.25, "lines"),
    legend.text          = element_text(margin = margin(r = 18)),
    plot.margin          = margin(t = 18, r = 10, b = 20, l = 10),
    panel.background     = element_rect(fill = "white", color = NA),
    plot.background      = element_rect(fill = "white", color = NA)
  )

# ---- Save --------------------------------------------------------------------
ggsave(out_file, p, width = 6, height = 4, dpi = 300, bg = "white")
message("✅ Wrote: ", normalizePath(out_file, winslash = "/"))

